using System;
using System.Collections.Generic;
using System.Globalization;
using System.IO;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Plotting;
using Atomic.Libraries.Plotting.Functions;

namespace Atomic.Thermodynamics.Electrons
{
	/// </summary>
	[Serializable]
	public class ElectronicStateDensity : IElectronicStateDensity
	{
		public ElectronicStateDensity(IEnumerable<ElectronicStateDensityPoint> points, double fermiEnergy, double electrons)
		{
			// Here the spin-polarized points are left undefined (null).
			Points = new ReadOnlyList<ElectronicStateDensityPoint>(points);
			FermiEnergy = fermiEnergy;
			Electrons = electrons;
		}

		public ElectronicStateDensity(IEnumerable<ElectronicStateDensitySpinPolarizedPoint> points, double fermiEnergy, double electrons)
			: this((IEnumerable<ElectronicStateDensityPoint>)points, fermiEnergy, electrons)
		{
			SpinPolarizedPoints = new ReadOnlyList<ElectronicStateDensitySpinPolarizedPoint>(points);
		}

		/*/// <summary>
		/// Parses VASP's DOSCAR file.
		/// </summary>
		public static VaspElectronicStateDensity Parse(double fermiEnergy, double electrons, TextReader reader)
		{
			List<VaspElectronicStateDensityPoint> points = new List<ElectronicStateDensityPoint>();

			// Skip header lines.
			for (int i = 0; i < 5; i++)
			{
				reader.ReadLine();
			}

			// Except last header line containing higher-accuracy values for energies.
			string line0 = reader.ReadLine().Trim();
			while (line0.Contains("  "))
			{
				line0 = line0.Replace("  ", " ");
			}

			string[] cols0 = line0.Split(' ');

			double emax = double.Parse(cols0[0], CultureInfo.InvariantCulture);
			double emin = double.Parse(cols0[1], CultureInfo.InvariantCulture);
			double n = int.Parse(cols0[2]);

			int j = 0;
			string line;
			while ((line = reader.ReadLine()) != null)
			{
				line = line.Trim();
				while (line.Contains("  "))
				{
					line = line.Replace("  ", " ");
				}

				string[] cols = line.Split(' ');

				if (cols.Length == 3)
				{
					//double e = double.Parse(cols[0], CultureInfo.InvariantCulture);
					double e = emin + (emax - emin) / (n - 1) * j;
					double d = double.Parse(cols[1], CultureInfo.InvariantCulture);
					double s = double.Parse(cols[2], CultureInfo.InvariantCulture);
					points.Add(new VaspElectronicStateDensityPoint(e, d, s));
				}
				else if (cols.Length == 5)
				{
					//double e = double.Parse(cols[0], CultureInfo.InvariantCulture);
					double e = emin + (emax - emin) / (n - 1) * j;
					double d1 = double.Parse(cols[1], CultureInfo.InvariantCulture);
					double d2 = double.Parse(cols[2], CultureInfo.InvariantCulture);
					double s1 = double.Parse(cols[3], CultureInfo.InvariantCulture);
					double s2 = double.Parse(cols[4], CultureInfo.InvariantCulture);
					points.Add(new VaspElectronicStateDensitySpinPolarizedPoint(e, d1, d2, s1, s2));
				}
				else
				{
					throw new NotSupportedException();
				}

				j++;
			}

			// First entry is used to determine whether the file contains spin-polarized information.
			if (points.Count > 0 && points[0] is VaspElectronicStateDensitySpinPolarizedPoint)
			{
				return new VaspElectronicStateDensity(points.Select(p => (VaspElectronicStateDensitySpinPolarizedPoint)p), fermiEnergy, electrons);
			}
			else
			{
				return new VaspElectronicStateDensity(points, fermiEnergy, electrons);
			}
		}*/

		public void Plot()
		{
			if (SpinPolarizedPoints != null)
			{
				IPlotFunction d1 = new LinearInterpolationFunction(SpinPolarizedPoints.Select(p => new PlaneVector(p.Energy, p.SpinUpDensity)));
				IPlotFunction d2 = new LinearInterpolationFunction(SpinPolarizedPoints.Select(p => new PlaneVector(p.Energy, p.SpinDownDensity)));
				new Gnuplot()
					.SetXLabel("Energy per unit cell, eV")
					.SetYLabel("Density, eV^{-1}")
					.Add(string.Format(CultureInfo.InvariantCulture, "set arrow from {0},0 to {0},{1} nohead", FermiEnergy, SpinPolarizedPoints.Max(p => Math.Max(p.SpinUpDensity, p.SpinDownDensity))))
					.Plot(
						new FunctionPlot(d1, SpinPolarizedPoints.Min(p => p.Energy), SpinPolarizedPoints.Max(p => p.Energy), LineColor.Green, new PlotLabel("d_up")),
						new FunctionPlot(d2, SpinPolarizedPoints.Min(p => p.Energy), SpinPolarizedPoints.Max(p => p.Energy), LineColor.Red, new PlotLabel("d_down"))
					)
					.Pause()
					.Run();
			}
			else
			{
				IPlotFunction d = new LinearInterpolationFunction(Points.Select(p => new PlaneVector(p.Energy, p.Density)));
				new Gnuplot()
					.SetXLabel("Energy per unit cell, eV")
					.SetYLabel("Density, eV^{-1}")
					.Add(string.Format(CultureInfo.InvariantCulture, "set arrow from {0},0 to {0},{1} nohead", FermiEnergy, Points.Max(p => p.Density)))
					.Plot(new FunctionPlot(d, Points.Min(p => p.Energy), Points.Max(p => p.Energy), LineColor.Black, new PlotLabel("d")))
					.Pause()
					.Run();
			}
		}

		/// <summary>
		/// Electronic density of states.
		/// </summary>
		public ReadOnlyList<ElectronicStateDensityPoint> Points
		{
			get;
			private set;
		}

		/// <summary>
		/// Spin-polarized electronic density of states.
		/// </summary>
		public ReadOnlyList<ElectronicStateDensitySpinPolarizedPoint> SpinPolarizedPoints
		{
			get;
			private set;
		}

		IEnumerable<IElectronicStateDensityPoint> IElectronicStateDensity.Points
		{
			get
			{
				return Points;
			}
		}

		public double FermiEnergy
		{
			get;
			private set;
		}

		public double Electrons
		{
			get;
			private set;
		}
	}
}
